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ABSTRACT 

We investigate the effects of neutrino heating and a-particle recombination on the hydrodynamics of core- 
collapse supernovae. Our focus is on the non-linear dynamics of the shock wave that forms in the collapse, 
and the assembly of positive energy material below it. To this end, we perform time-dependent hydrodynamic 
simulations with FLASH2.5 in spherical and axial symmetry. These generalize our previous calculations by 
allowing for bulk neutrino heating and for nuclear statistical equilibrium between n, p and a. The heating rate 
is freely tunable, as is the starting radius of the shock relative to the recombination radius of a-particles. An ex- 
plosion in spherical symmetry involves the excitation of an overstable mode, which may be viewed as the I = 
version of the 'Standing Accretion Shock Instability' . In two-dimensional simulations, non-spherical deforma- 
tions of the shock are driven by plumes of material with positive Bernoulli parameter, which are concentrated 
well outside the zone of strong neutrino heating. The non-spherical modes of the shock reach a large amplitude 
only when the heating rate is also high enough to excite convection below the shock. The critical heating rate 
that causes an explosion depends sensitively on the initial position of the shock relative to the recombination 
radius. Weaker heating is required to drive an explosion in two dimensions than in one, but the difference also 
depends on the size of the shock. Forcing the infalling heavy nuclei to break up into n and p below the shock 
only causes a slight increase in the critical heating rate, except when the shock starts out at a large radius. This 
shows that heating by neutrinos (or some other mechanism) must play a significant role in pushing the shock 
far enough out that recombination heating takes over. 

Subject headings: hydrodynamics — instabilities — nuclear reactions, nucleosynthesis, abundances — shock 
waves — supernovae: general 
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1. INTRODUCTION 

Although tremendous progress has been made on the mech- 
anism of core-collapse supernovae in recent years, we still 
do not have a clear picture of a robust path to an explosion 
in stars that form iron cores. Heating by the absorption of 
electron-type neutrinos significantly modifies the settling flow 
below the bounce sho ck but - in spite of early positive results 
dBethe & W ilson 1985) - explosions are obtained in spherical 
collapse calculati ons only if the pro genitor star is lighter than 
about 10-12 Mq dKitaura et al . 2006). M ore massive stars fa il 
to explode in spherical symmetry dLiebendo rfer et alj|2001l) . 

Two-dimensional collapse calculations show strong de- 
for mations of the shock and convective mo t ions below 
it dBurrows et all 11995b Uanka & Mueller! 1 1996t iBuras et alj 
l2006allbl:lBurrows et alj2006L 120071: iMarek & Jankall2009l) . It 
has long been recognized that convection increases the res- 
idency time o f settling material in the zone of strong neu- 
trino heating dHerant et alJl!992h . It is also becoming clear 
that multidimensional explosions require the assembly of a 
smaller amount of material with positive energy, but the de- 
tails of how this happens remain murky. 

An early treatment of shock breakout by iBethd d 19971) 
focused on the strong gradient in the ram pressure of the 
infalling material, but implicitly assumed that the shocked 
material had already gained positive energy. If large-scale 
density inhomogeneities are present below the shock, they 
will trigger a finite-amplitude, dipolar instability, thereby al- 
lowing accretion to continue simultaneously with the expan- 
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sion of positive-energy fluid (Thompson 2000). The ac- 
cretion shock is also capable of a dipolar oscillation which 
leads, above a critical amplitude, to a bifurcation between 
freshly infalling mate rial and material shocked at earlier times 
(Blondin et al. 2003). Such an oscillation is easily excited in 
a spherical flow composed of a zero-energy, poly tropic fluid 
(Blondin & Mezzacappa 2006) via a linear feedback between 
ingoing vortex and en tropy waves and outgoing sound waves 
(Fogli zzo et al.l 120071) . It can also be excited indirectly by 
neutrino heating, which if sufficiently str ong will drive large- 
scale buoyant motio ns below the shock dHerant et al J 1 1 994c 
Foglizzo et al. 2006). For relatively weak heating, the dipo- 
lar oscillation can decrease the damping effect of advection 
and trig ger convective motions that would other wise be sup- 
pressed (Foglizzo et al. 2006; Scheck et al. 2008). 

A significant sink of thermal energy in the accretion flow 
arises from the dissociation of heavy nuclei. The heavy el- 
ements that flow through the shock are broken up into a- 
particles and nucleons when exposed to the high temperature 
(> 1 MeV) of the postshock region. The Bernoulli parame- 
ter b of the shocked fluid then becomes substantially negative. 
A significant fraction of this dissociation en ergy can be r ecov- 
ered if nucleons recombine into a-particles (Bethe 1996]). But 
for this to occur, a decrease in the temperature is required and 
thus the shock must expand significantly beyond the radius at 
which it typically stalls (— 100- 150km). 

O ne of the primary goals here , and in a previous pa- 
per dFernandez & Th ompson 2009) [hereafter Paper I], is to 
gauge the relative importance of these effects in setting the 
stage for a successful explosion. The persistence and ampli- 
tude of a dipolar oscillation can only be reliably measured in 
fully three-dimensional simulations (there are preliminary in- 
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dications that it is less pr ominent in three spa tial dimensions 
than in axial symmetry; Iwa kami et al.l l2008). On the other 
hand, the interplay between a-particle recombination and hy- 
drodynamical instabilities has received little attention. Al- 
though recombination is certainly present in previous numeri- 
cal studies which employ finite-temperature equations of state 
(EOSs), it should be kept in mind that considerable uncertain- 
ties in the EOS remain at supranuclear densities. A softening 
or hardening of the EOS feeds back on the position of the 
shock for a given pre-collapse stellar model ( Ma rek & Jankal 
2009). The parametric study of th e critical neutrino luminos- 
ity by Murphy & Burrows (2008) is based on a single EOS; 
they obtain explosions in which the shock seems to break out 
from nearly the same radial position at ~ 250-300 km. Vari- 
ations in the density profile of the progenitor star will simi- 
larly modify the position of the shock, the concentration of 
a-particles below it, and the critical neutrino luminosity for 
an explosion. 

In this paper, we study the interplay between non-spherical 
shock oscillations, neutrino heating, and a-particle recombi- 
nation, when the heating rate is pushed high enough for an 
explosion to occur. Our focus is on the stalled shock phase, 
between ~ 100 ms and 1 s after bounce. In a one-dimensional 
calculation, the accretion flow reaches a quasi-steady state 
during this interval, and the shock gradua lly recedes (e.g. 
iLiebendorfer et alJl200ltlBuras et alj|2006ah . 

Our approach is to introduce EOSs of increasing complex- 
ity into one- and two-dimensional, time-dependent hydrody- 
namic simulations. To this end, we use the code FLASH2.5 
dFryxell et a l. 2000), which is well tested in probl ems involv- 
ing nu clear energy release in compressible fluids ( Cal der et"aT] 
2002). We adopt a steady state model as our initial condition, 
and a constant mass accretion rate, neutrino luminosity, and 
fixed inner boundary. The steady-state appr oximation to the 
stalled shock phase was first introduced by Burrows & Goshy 
(119931) . and has recently been used by Ohnis hiet alJ (12006) to 
study the non-linear development of the shocked flow with a 
semi-realistic equation of state and neutrino heating. In Paper 
I we modeled the accretion flow as a polytropic fluid, from 
which a fixed dissociation energy is removed immediately be- 
low the shock, and allowed for neutrino cooling but not heat- 
ing. 

Here we generalize this model to allow both for heating, 
and for nuclear statistical equilibrium (NSE) between neu- 
trons, protons, and a-particles. Heating is introduced in a 
simple, parametrized way, without any attempt at simulating 
neutrino transport. The nuclear abundances are calculated as 
a function of pressure p and density p, using a complete finite- 
temperature, partially degenerate EOS. Our model for the 
shocked material retains one significant simplification from 
Paper I: we do not allow the electron fraction Y e to vary with 
position below the accretion shock. Here there are two com- 
peting effects: electron captures tend to reduce Y e , whereas 
absorption of v e and v e tends to drive high-entropy material 
below the shock toward Y e ~ 0.5. Since we are interested es- 
pecially in the dynamics of this high-entropy material, we set 
Y e = 0.5 when evaluating the a-particle abundance. To obtain 
a realistic density profile, we continue to approximate the in- 
ternal energy of the fluid as that of a polytropic fluid with a 
fixed index 7 = i. In reality, the equation of state between the 
neutrinosphere and the shock depends in a complicated way 
on the degeneracy of the electrons and the effects of electron 
captures. The consequences of introducing these additional 



degrees of freedom will be examined in future work. 

In spite of these simplifications, our results already show 
many similarities with more elaborate collapse calculations. 
Spherical explosions are due to a global instability resem- 
bling the one-dimensional Standing Accretion Shock Insta- 
bility (SASI), but modified by heating. As in Paper I, we 
find that the period of the £ = mode remains close to twice 
the post-shock advection time. Strong deformations of the 
shock in two-dimensional runs are driven by material with 
positive Bernoulli parameter, which generally resides outside 
the radius r a where the gravitational binding energy of an a- 
particle is equal to its nuclear binding energy. The recombina- 
tion of a-particles plays a major role in creating this positive- 
energy material, but for this to happen the shock must be 
pushed beyond ~ 200 km from the neutronized core. 

In this paper, we consider only neutrino heating as the im- 
petus for the initial expansion of the shock, rather than more 
exotic effects such as rotation or magnetic fields. We find that 
the critical heating rate is a strong function of the initial po- 
sition of the shock with respect to r a , which implies that a 
much higher neutrino luminosity is needed to revive a shock 
that stalls well inside r a . The difference in the critical heating 
rate between one- and two-dimensional simulations also de- 
pends on the size of the shock, and thence on the structure of 
the forming neutron star. 

We find that the shock develops a dipolar oscillation with a 
large amplitude only when the heating rate is also high enough 
to trigger a strong convective instability. We therefore sur- 
mise that buoyant motions driven by neutrino heating play a 
major role in driving the dipolar oscillations that are seen in 
more complete simulations of core collapse. Some evidence 
is found that acoustic wave emission by the convective mo- 
tions also can play a role. We investigate the possibility of a 
heat engine within the gain layer (the region with a net excess 
of neutrino heating over cooling). We find that most of the 
heat deposition by neutrinos is concentrated in lateral flows at 
the base of the prominent convective cells. At the threshold 
for an explosion, neutrino heating plays a key role in pushing 
material to positive Bernoulli parameter, but only if the shock 
starts well inside r a . 

The plan of the paper is as follows. Section |2]presents our 
numerical setup, treatment of heating, cooling, and nuclear 
dissociation, and outlines the sequences of models. Sections 
[3] and |4] show results from one- and two-dimensional simula- 
tions, respectively. We focus on the relative effectiveness of 
neutrino heating and a-particle recombination in driving an 
explosion, and the relation between Bernoulli parameter and 
large-scale deformations of the shock. The critical heating 
rate for an explosion is analyzed in §|5] and the competition 
between advective-acoustic feedback and convective instabil- 
ity is discussed in ^6] We summarize our findings in §|7] The 
appendices contain details about our EOS and the numerical 
setup. 

2. NUMERICAL MODEL 

As in Paper I, the initial configuration is a steady, spheri- 
cally symmetric flow onto a gravitating point mass M. The 
flow contains a standing shock wave, and the settling flow 
below the shock cools radiatively in a narrow layer outside 
the inner boundary of the simulation volume. The space of 
such models is labeled basically by three parameters: accre- 
tion rate M, luminosity L v in electron neutrinos and anti- 
neutrinos, and the radius r* of the base of the settling flow, 
which corresponds roughly to the neutrinosphere radius. The 
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mass M of the collapsed material represents a fourth param- 
eter, but it covers a narrower range than the other three. The 
infalling material is significantly de-leptonized before hitting 
the shock only in the fir st 50 ms or so of the collapse (e.g. 
lLiebendorfefetal]|200l . 

We explore a two-dimensional surface through this three- 
dimensional parameter space by i) fixing the ratio of r* to the 
initial shock radius r,o in the absence of heating (r»/r,o = 0.4); 
ii) allowing r s o to vary with respect to an appropriately chosen 
physical radius; and then iii) increasing the level of heating 
until an explosion is uncovered. In the full problem, the shock 
radius at zero heating is a unique function of M and r*, with 
a small additional dependen ce on M and the composi tion of 
the flow outside the shock dHouck & Chevalier! [l 9921) . The 
secular cooling of the collapsed core forces a gradual decrease 
in r*, and M also varies with time and with progenitor model. 

Given the important role that a-particle recombination 
plays in the final stages of an explosion, we implement ii) by 
referencing r& to the radius where the gravitational binding 
energy of an a-particle equals its nuclear binding energy, 



GMm C) 
Q a 



: 254Mi.3 km. 



(1) 



Here Q a ~ 28.30 MeV is the energy needed to break up an 
a-particle into 2n and 2p, m a the mass of an a-particle, and 
Mi .3 = M /(\3Mq). Choice i) allows us to consider models 
that have, implicitly, both a range of physical values of and 
a range of M. It is, of course, made partly for computational 
simplicity (the limited size of the computational domain) and 
also to facilitate a comparison between models that have dif- 
ferent values of r s o/r a . Nuclear dissociation is taken into ac- 
count either by removing a fixed specific energy s right below 
the shock, or by enforcing NSE between n, p, and a through- 
out the settling flow. Once this choice is made, the normal- 
ization of the cooling function is adjusted to give r»/r s o = 0.4. 
The heating rate remains freely adjustable thereafter. 

We adopt this simplification because we do not intend to 
find the precise value of the critical neutrino luminosity, but 
instead to probe the behavior of the system around this critical 
point, whatever its absolute value. 

We now describe the key components of this model in more 
detail, and explain the setup of the hydrodynamic calcula- 
tions. As in Paper I, the time evolution is carried out using the 
second-ord er, Godunov-type, a daptive-mesh-refinement code 



secona-oraer, Uoaunov-type, aa 
FLASH2.5 (IFrvxell et al.11 2000). 



2.1. Initial Conditions 

The introduction of heating causes a change in the structure 
of the initial flow configuration. The radius r s of the shock in 
the time-independent solution to the flow equations increases 
with heating rate; that is, r s > r s o- The material above the 
shock is weakly bound to the protoneutron star, and in practice 
can be taken to have a zero Bernoulli parameter 



1 - 7 P GM 
b = —v H . 

2 7— 1 p r 



(2) 



Here v is the total fluid velocity, which is radial in the initial 
condition, p is the pressure, p is the mass density, and G is 
Newton's constant. The flow upstream of the shock is adia- 
batic and has a Mach number M. \ = 5 at a radius r = r s o. 

The composition of the fluid is very different upstream and 
downstream of the shock. Changes in internal energy due to 
nuclear dissociation and recombination are taken into account 



using the two models described in 92.1.11 For the internal 
energy density of the fluid <?, we continue to use the polytropic 
relation e = p/i'j-l); hence the second term on the right-hand 
side of equation (fJJ. Because we are not explicitly including 
changes in electron fraction due to weak processes, we keep 
7 = 4/3 for both models of nuclear dissociation. This largely 
determines the density profile inside a radius ~ \r a , where 
the a-particle abundance is very low. 

The upstream and downstream flow profiles are connected 
through the Rankine-Hugoniot jump conditions, which are 
modified so as to allow a decrement e in b across the shock. 
The resulting compression factor is (Paper I) 



K =^ = (7+l) 
Pi 



(7+ -MI 2 ) 



fl-M7 2 ) 2 + (7 2 -D^ 



(3) 



which reduces to k — > (7+ l)/(7~ 1) for Mi — » 00 and e = 
0. Throughout this paper, the specific nuclear dissociation 
energy e is defined to be positive. The subscripts 1 and 2 
denote upstream and downstream variables, respectively. 

All flow variables are made dimensionless by scaling radii 
to r s o, velocities to Vffo = (2GM / r S Q) x l 2 , timescales to fffo = 
'so/vffo, a nd de nsities to the upstream density at r= r s o, pi(r s o) 
[equation IB 131 . See Paper I for further details. Throughout 
the paper we denote the average of a function F(X, ...) over 
some variable X by (F) X - 

2.1.1. Nuclear Dissociation 

We model nuclear dissociation in two ways. First, we re- 
move a fixed specific energy e right below the shock, as done 
in Paper I. This represents the prompt and complete breakup 
of whatever heavy nuclei are present in the upstream flow. 
The main limitation of this approximation is that the dissoci- 
ation energy does not change with the radius (or inclination) 
of the shock. The main advantage is simplicity: e is indepen- 
dent of any dimensional parameters and can be expressed as a 
fraction of v^ . 

We also use a more accurate dissociation model which 
allows for NSE between a-particles and nucleons. 3 Dur- 
ing the stalled shock phase of core-collapse supernovae, the 
shock sits at r ~ 100-200 km, with a postshock temperature 
T > 1 MeV and density p > 10 9 g cm" 3 . In these conditions, 
the heavy nuclei flowing through the shock are broken up into 
a, p, and n. 

A range of isotopes are present in the iron core of a mas- 
sive s tar as well as in nuclear burning shells dWooslev et alj 
2002), but since the binding energy per nucleon varies only 
by ~ 10% we simply assume a single type of nucleus in the 
upstream flow. We focus here on the later stages of the stalled 
shock phase, during which the oxygen shell is accreted. An 
energy Qq ~ 14.44 MeV must be injected to dissociate an 
16 nucleus into 4 a-particles dAudi et a l. 2003), which cor- 
responds to the specific dissociation energy 



£o = ^0.038^3 
mo 



150 km 



y|(r). 



(4) 



3 Although heavier nuclei can begin to recombine once the shock moves 
significantly beyond r a , this generally occurs only after the threshold for an 
explosion has been reached, and makes a modest additional contribution to 
the recombination energy. 
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Here mo ~ 16m„ is the mass of an oxygen nucleus, with m u 
the atomic mass unit. The smallness of this number indicates 
that little oxygen survives in the post-shock flow, and so we 
set the equilibrium mass fraction of oxygen to zero below 
the shock, X Q q = 0. The binding energy of an a-particle is 
of course much larger, giving 



: 0.297M7. 3 



150 km 



VffO). 



(5) 



We find that a-particles appear in significant numbers only 
at relatively large radii (> 0.5 r a ) and in material that has ei- 
ther i) been significantly heated by electron neutrinos closer 
to the neutrinosphere; or ii) been freshly shocked outside r a . 
The electrons are only mildly degenerate in material that has 
a high entropy and a-particle content, so that neutrino heat- 
ing drives Y e close to ~ 0.5 (or even slightly above: see, e.g., 
IBuras et al]|2006bl) . We therefore set Y e = 0.5 in the Saha 
equation that determines the equilibrium mass fractions X^ q , 
Xp q and X^ = l-X^ q -Xp q . These quantities are tabulated 
as functions of p and p using an ideal, finite-temperature and 
partially degenerate equation of state for electrons and nucle- 
ons; see Appendix lAl for details. Specific choices must then 
be made for the parameters r s o, M, and M; we generally take 
M= 1.3M Q andM = 0.3M© s" 1 , but allow r s0 to vary. An in- 
vestigation of how changes in Y e feed back onto the formation 
of a-particles is left for future work. 

A specific energy 



e n uc = -X (e + e a )- (X a -X sq [p,p\) e D 



(6) 



is either released or absorbed within a single time step (it can 
be of either sign). Here Xo is non-vanishing only for fluid ele- 
ments that have just passed across the shock, and we have set 
X^ = 0. The quantity © is introduced as an energy source 
term in FLASH, and from it one readily obtains a rate of re- 
lease of nuclear binding energy per unit mass, 



At 



At ' 



(7) 



where Af is the simulation time step. 

In the initial condition, the dissociation energy at the shock 
is obtained from equation (|6]l using Xo = 1 and X a = up- 
stream of the shock: 



e(t = Q) = eo+(l-X^ q [p 2 ,P2])e 



(8) 



FigureQ]shows how e(f = 0) and X^ depend on the shock ra- 
dius r s o, for upstream flows composed 4 of pure ie O and 56 Fe, 
and for different values of M. The dissociation energy is ap- 
proximately constant inside ~ 75 km, where the downstream 
flow is composed of free nucleons, but decreases at greater 
distances, remaining ~ 40% of the gravitational binding en- 
ergy at the shock. The mass fraction of a-particles reaches 
50% at r = 150-175 km, with a weak dependence on M. 

2.1.2. Heating and Cooling in the Post-shock Flow 

To allow direct comparison with our previous results, we 
employ a cooling rate per unit volume of the form 



se c = c P a p h 



(9) 



In the case where the upstream flow is pure Fe, we replace eo in equa- 



tion {8) with ep ( 
tron fraction to Y, 



■■ Qfe/mpe ~ 0.093^3 (r/150 km)v|(r), and set the elec- 
= 26/56 in the NSE calculation behind the shock. 
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FIG. 1. — Equilibrium mass fraction of a-particles X„ , and ratio of ini- 
tial dissociation energy e(t = 0) [equation [H] to vjL behind a spherical shock 
positioned at radius r s o. Curves of different shadings correspond to dif- 
ferent mass accretion rates. The Rankine-Hugoniot shock jump conditions 
and dissociation energy are calculated self-consistently, as described in Ap- 
pendix|E] Square brackets refer to the upstream composition of the accretion 
flow, which for simplicity is taken to be pure 56 Fe or 16 0. The Mach number 
upstream of the shock is M\ = 5, and the central mass is M = 1.3M0. 



with a = 1.5, b = 2.5, and C a normalization constant. As in 
Paper I, we include a gaussian entropy cutoff to prevent run- 
away cooling. The exponents in equation (0 represent cool- 
ing dominated by the capture of relativistic, n on-degenerat e 
electrons and positrons on free nucleons (e.g.. lBethall990h . 
Inside the radius where the electrons become strongly de- 
generate, and a-particles are largely absent, one has Jz?c oc 

3/2 3 

p e n p oc (Y e p) . This gives essentially the same dependence 
of .5fc on r as equation © when Y e = constant and 7=3 (cor- 
responding to p oc r~ 3 in a nearly adiabatic settling flow). In 
more realistic collapse calculations, Y e grows with radius be- 
tween the neutrinosphere and the shock, but p tends to de- 
crease more rapidly than ~ r" 3 (e.g. IBuras et al] |2006a). Our 
chosen form for the cooling function results in a slightly wider 
gain region and, therefore, a slightly lower critical heating rate 
for an explosion. The bulk of the cooling occurs in a narrow 
layer close to the accretor at r = r*, and the accreted material 
accumulates in the first few computational cells adjacent to 
the inner boundary without a major effect on the rest of the 
flow. 

We model neutrino heating as a local energy generation rate 
per unit volume of the form 



J? H = H(l-X a )p/r 2 



(10) 



The normalization constant H measures the strength of the 
heating. The factor (1 —X a ) accounts for the fact that the cross 
section for neutrino absorpti on by a-part icles is much smaller 
than that for free nucleons (lBethdll990i) . For simplicity, we 
do not include the flux factor due to the transition between 
diffusion and free-streaming. Our focus here is on the nature 
of the instabilities occurring in the flow near the threshold for 
an explosion, and we do not attempt a numerical evaluation of 
the critical heating rate. 

An additional energy source term arises from the change in 
the equilibrium fraction of a-particles as they are advected in 
the steady state initial solution. The instantaneous adjustment 
of X a to its equilibrium value, combined with equation (O, 
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yields an energy generation rate per unit volume 
dX? 



J£ a = pVE a - 



dr 



: pve a 



dp dX? dp 



dp dr dp dr 



(11) 



This energy generation rate is negative, as the temperature 
increases inwards and thus the a-particle fraction decreases 
with decreasing radius (v is negative). 

2.1.3. Numerical Setup 

In our time dependent calculations, we use one-dimensional 
and two-dimensional spherical coordinates with baseline res- 
olution Arb aS e = r s o/320 and Af?b ase = 7r/192, with one extra 
level of mesh refinement inside r = r* + 0.1(r s o-r*) to better 
resolve the steep density gradient that arises in the cooling 
layer. We do not employ a hybrid Riemann solver because we 
do not see th e appearance of the odd-even decoupling insta- 
bility (Ou ir^[l99l . 

We employ a reflecting inner boundary condition at r = r* 
for the sake of simplicity; we do not attempt to model the 
protoneutron star (as done by Murphv & Burrows 2008) or 
its contraction throu gh a mo ving inner bound ary [as done 
bv lSchecketal] d2006l) and iScheck etail ( 120081) 1. The outer 
boundary condition is kept fixed at r = 7r s o, and is set by the 
upstream flow at that position. 

To trigger convection below the shock, we introduce ran- 
dom cell-to-cell velocity perturbations in v, and v$ at t = 0, 
with an amplitude 1 % of the steady state radial velocity. To 
study the interplay between shock oscillations and convection, 
we also drop overdense shells with a given Legendre index I, 
as done in Paper I, without random velocity perturbations. 

In order to track the residency time of the fluid in the gain 
region, we assign a scalar to each spherically symmetric mass 
shell in the upstream flow. This scalar is passively advected 
by FLASH2.5. Through this technique, we are able to assign 
a "fluid" time to each element in the domain, corresponding 
to the time at which the mass shell would cross the instan- 
taneous angle averaged shock position if advected from the 
outer boundary at the upstream velocity: 



f 0B dr 
t F = t B+ I i— r- 



(r,(t))e H 



(12) 



Here ?ob is the time at which the fluid enters through the outer 
radial boundary at r = tob, and (r s (f))e is the angle averaged 
shock position. Initially, ?ob = and all the fluid below the 
shock is set to tf = 0. This prescription works well for statis- 
tical studies ( 134.21 ). tracing large scale fluid patches, despite 
some inevitable turbulent mixing on small scales. 

An explosion is defined as either i) a collision between 
the shock and the outer boundary of the simulation volume 
(r = 7r s o) within lOOOffto of the start of the simulation; or ii) 
in the special case of the one-dimensional constant-e mod- 
els, a transient expansion that breaks a quasi-steady pattern 
within the same timeframe. Even in the spherically symmet- 
ric simulations, very small changes in heating rate can lead to 
dramatic changes in shock behavior, and so this definition of 
explosion is good enough for our purposes. 

2.2. Model Sequences 

We choose six sequences of models, each with a range of 
heating parameters H > 0, and each evolved both in spher- 
ical and axial symmetry. Their parameters are summarized 
in Table Q] In each sequence, the normalization of the cool- 
ing function is chosen so that r*/r s = 0.4 at zero heating. 




r/r s o 

FIG. 2. — Sample initial density profiles, which are solutions to the spher- 
ically symmetric and time-independent flow equations. Panel (a) shows the 
zero-heating configurations for all the sequences listed in Table [Q Other 
parameters are {7 = 4/3, M.\(r s ia) = 5} for all configurations, and {M = 
O.3M s -1 , M= 1.3M } for the NSE models. Panel (b) shows a sequence 
with a fixed cooling function and range of heating rates (H is given in units 
of r s oVj fo ). The dashed line shows the upstream flow. Panel (c) shows a 
sequence with different equations of state, r,o = 125 km, and H = 0. The la- 
bels "+o" and "—a" mean with and without a-particles included in the EOS, 
while "full" means that the EOS explicitly includes finite-temperature and 
partially degenerate electrons, black body photons, and ideal-gas ions. All 
other parameters are the same as in (a). Only the 7 = 4/3 upstream flow is 
shown. See Appendix IBlfor further details. 



Three sequences have a constant dissociation energy, which 
take the values e/vg = {0.1,0.15,0.2}. The other three se- 
quences assume NSE below the shock, and have shock radii 
r s o = {50,75, 125} km at zero heating. This means that the 
physical value of the cooling radius also takes on different val- 
ues, namely {20,30,50} km. In effect, our models are prob- 
ing different sizes for the neutrinosphere, and different times 
following the collapse. The other parameters in the NSE mod- 
els are M = 1.3M Q and M = 0.3M Q s" 1 . 

Table[T]samples some properties of a few models from each 
sequence: one with zero heating, another with H close to the 
critical value for an explosion, and a third with the largest 
heating parameter that will allow a steady solution. Note that 
the shock starts out at ~ 1 .3 r s o in the time-independent, spher- 
ical flow solution, and quickly saturates at ~ (1.8— 2)r s o in the 
two-dimensional models with heating just below threshold for 
an explosion. The quantity ejv\ references the dissociation 
energy to (twice) the kinetic energy of the upstream flow, and 
is the key free parameter determining the compression rate k 
across the shock (equation Q). 

When examining how the prescription for nuclear dissoci- 
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ation influences the results, we will focus on the e = 0.15vj? t0 
sequence and the NSE sequence with r su = 75 km, which have 
similar initial density profiles (due to the low initial a-particle 
abundance in the NSE model). 

The six initial models at zero heating are shown in Fig- 
ure 12^. Panel (b) shows the sequence of initial models with 
r S Q = 75km and a range of heating parameters. The model 
with H = 0.007vj? f0 r s o is close to the threshold for an explo- 
sion, while the one with H = 0.009v^ r su is well above thresh- 
old. At higher values of H, cooling by a-particle dissociation 
(equation fTTI ) can be significant in a layer below the shock, 
causing the density profile to steepen slightly. 

Fig [5J; shows how our constant-7, ideal gas approximation 
to the internal energy of the flow compares with the full EOS 
containing finite-temperature and partially degenerate elec- 
trons (see Appendix IB1 for details). The curves labeled 
include our prescription for heating/cooling by a-particle re- 
combination/dissociation, and those labeled "-a" do not. We 
show the sequence with the largest shock radius (r s o =125 
km) so that NSE allows some as to be present. The neglect 
of electron captures below the shock results in an adiabatic 
index between 4/3 and 5/3 in the zone where a-particles are 
absent. This causes the EOS to stiffen, so that the density pro- 
file is well approximated by an ideal gas with 7 ~ 1 .48 at zero 
heating. Adding in heating tends to flatten the density profile 
even more, and with 7 = 1 .48 it would be much flatter than 
is typically seen in a realistic core collapse model. Hence we 
choose an EOS with 7 = 4/3. 

3. ONE DIMENSIONAL SIMULATIONS 

3.1. Shock Oscillations and Transition to Explosion 

An explosion in spherical symmetry involves the develop- 
ment of an unstable £ = SASI mode. We showed in Paper 
I that, in the absence of neutrino heating, the period of this 
mode is essentially twice the post-shock advection time. As 
heating is introduced into the flow, we find that this relation is 
maintained. The £ = mode is damped until the heating rate 




0.004 



0.01 



FIG. 3. — Linear growth rates (top) and oscillation frequencies (bottom) of 
one-dimensional models with constant dissociation energy e, as a function of 
heating parameter H around the threshold for explosion. Stars denote con- 
figurations that explode within 1000%o- Increased heating makes the system 
more unstable because the density profile flattens, akin to an increase in 7. 
Dotted lines show the frequency u) ox = 27r/(2f a d v ). Oscillation frequencies 
decrease with increasing heating rate beca use r s moves out relative to r„ , so 
that the advection time f aI j v (equation |13| ) increases. Increasing the dissoci- 
ation energy raises the oscillation period, and so a somewhat higher heating 
rate is required to obtain an explosion in a finite interval. 



is pushed above a critical value, which we now discuss. 

It should be emphasized that thi s critical heating rate is gen - 
erally lower than that defined by iBurrows & Goshvl (11993b . 
which marked the disappearance of a steady, spherically sym- 
metric solution to the flow equations. Large amplitude shock 
oscillations in spherical symmetry have been witnessed near 
the thr esh old for explosion in calc ulations by Ohnis hTet alj 
(2006) and Murphy & Bu rrows! (120081) . Both calculations em- 
ployed a realistic EOS, but like us included neutrino heating 
as a local source term in the energy eq uation. Oscillations 
have also been seen bv lBuras et alj d2006bl) in more elaborate 
calculations with Boltzmann neutrino transport. 

The origin of the spherically symmetric SASI oscillation 
can be briefly summarized as follows. An initial outward 
shock displacement generates an entropy perturbation, which 
is negative for 7 < 5/3. This entropy perturbation is advected 
down to the cooling layer, where it causes, at constant am- 
bient pressure, an increase in the cooling rate, SJ^cl-^c = 
— [(7 — l)/"f](b- a)5S > 0. The resulting negative pressure 
perturbation is rapidly communicated to the shock, which re- 
cedes and generates an entropy perturbation of the opposing 
sign. One more iteration results in a shock displacement of 
the same sign as the initial displacement, and allows the cycle 
to close. The duration of the £ = mode is therefore nearly 
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FIG. 4. — Shock radius as a function of time for two sequences of one- 
dimensional simulations. The upper panel shows runs with constant dissocia- 
tion energy e = 0.15vjLj, and a range of heating coefficients H near the critical 
value H a = (0.006625 ± 0.000125)vJ fo r s0 . The lower panel shows runs with 
r s o = 75 km that include recombination of a-particles. In this case, is ini- 
tially negligible everywhere below the shock (see Table[T}, but grows as the 
shock expands. The horizontal dotted line labels the radius r a at which the 
nuclear binding energy Q a of an a-particle equals its gravitational binding 
energy (equation [TJ. The critical heating for this second sequence is lower, 
Ha = (0.006125 ± 0.000125)v| r 8 o. 

twice the advection time from the shock to the cooling layer, 

-= 2 fn- < 13 > 

^osc J r , N 

The cycle is stable for 7 = 4/3 and r*/V s o = 0.4. 

When heating is added, the density profile flattens. Increas- 
ing 7 has the same effect, and has been found to push up 
the growth rate of linear SASI modes (e.g. Paper I). There 
is a critical heating rate for which the damping effect of the 
spherically symmetric SASI is neutralized and there is no net 
growth. We find that, once the heating rate exceeds this criti- 
cal value, the system always explodes. 

We therefore define the critical heating rate in our spheri- 
cally symmetric simulations to be the minimum heating rate 
for growing shock oscillations. 5 Figure [3] shows real and 
imaginary eigenfrequencies as a function of heating rate for 
our one-dimensional initial configurations with constant e. 
T he curves were obtain ed by solving the differential system 
of iFoglizzo et al.l (120071) . modified to account for a constant 
rate of nuclear dissociation (Paper I) as well as incorporating 

5 We define our critical heating parameter H a to be the average of the val- 
ues in the exploding and non-exploding runs that are closest to the threshold 
for explosion, within our fiducial 1000fffo cutoff. 
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FIG. 5. — Radial profiles of various quantities during shock breakout in 
the NSE model with H = \.(BH a and r s0 = 75 km (Figure 0. Top panel: 
mass fraction of a-particles. Second panel: rate of release of specific nu- 
clear binding energy de nuc /df compared with the (adiabatic) rate of change 
of enthalpy Wad f equation 1141 . Third panel: net neutrino heating rate per 
unit volume —J?c (thin solid curves) and de nuc /df (thick dashed curves), 
both normalized to the local value of c 2 s = yp/ p. Bottom panel: radial veloc- 
ity normalized to Vffo at radius r s Q. Both de nuc /d/ and w a( j are smoothed in 
radius for clarity. 

the heating function in equation ( TTOb . The runs marked by 
stars explode within a time 1000fffo, and so require a small, 
but finite, positive growth rate. 

In an exploding run, the expansions become longer and con- 
tractions shorter as the shock oscillation develops a large am- 
plitude. Eventually the accretion flow is halted during a con- 
traction. This marks the point of explosion, beyond which the 
feedback between the shock and the cooling layer is broken. 
Material then tends to pile up in the gain region, is further 
heated, and more material reverses direction. The net effect 
is to push the shock outward. Movie 1 in the online material 
illustrates this chain of events. 
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FIG. 6. — Left panel: Angle-averaged shock radius (solid lines) and maximum shock radius (dotted lines) for various two-dimensional models around the 
threshold for explosion. Upper panel shows runs with constant dissociation energy e = 0J_5vj f0 , while lower panels displays NSE runs with r s o as labeled. 
Critical heating rates H CI are different for each configuration, and can be found in Figure 1151 Right panel: Black lines show expansion timescale of maximum 
shock radius ? C xp ~ r s , ma x/|dr s , max /d/ 1, computed using a polynomial fit for the runs just above the threshold for explosion (corresponding to the black dashed 
lines on th e left panels). Green and red lines show the average residency time over the 50% and 10% of the gain region volume with highest ; rcs , respectively 
(see j|4.2| for the definition of this timescale). Shock breakout occurs whenever f cxp ~ (? re s)vol> except in the model where recombination heating is dominant 
(r s0 = 125 km). 
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3.2. Effects of Alpha-Particle Recombination 

Shock breakout is controlled by the build-up of positive en- 
ergy fluid downstream of the shock, and therefore is sensitive 
to the density profile below the shock. Heating by neutrinos 
is concentrated fairly close to the protoneutron star, inside a 
distance ~ (2-3)r». Heating by a-particle recombination is 
concentrated at a greater distance ~ r a (equation H]), but still 
can reach a comparable amplitude. 

The dependence of shock breakout on heating rate is dis- 
played in Figure|4]for two accretion models and several values 
of H close to H a (see Table [TJ. The initial expansion of the 
shock during the explosion phase is very similar for models 
with constant e and with NSE in the shocked fluid. However, 
the time evolution bifurcates near the radius r a . 

Figure [5] shows successive profiles of the shocked flow in 
the exploding run with H = 1 .09H a and r& = 75 km. The 
a-particle fraction approaches unity as the shock reaches the 
radius r a . The second panel shows the specific nuclear energy 
generation rate [equation (O] normalized to the adiabatic rate 
of change of the enthalpy, 

w ad = -^=-c?V-v. (14) 
p at 

Here c s = (jp/p) 1 ^ 2 is the sound speed. The third panel com- 
pares the amplitude and distribution of neutrino and recombi- 
nation heating, and the bottom panel plots the radial velocity 
in the postshock region. 

We can summarize this behavior as follows: during the ini- 
tial expansion phase, fluid below the shock continues to move 
inward, and the dissociation of a-particles removes energy 
from the flow (as expected from equation ifTTI ). Some fluid 
behind the shock begins to move outward around 300%), but 
nuclear dissociation still causes a net loss of internal energy. 
However, the recombination of a-particles sets in above r a , 
especially in regions where X a < 0.5. By the time the shock 
hits the outer boundary, de nuc /df exceeds one-half of |w a d|- 

The dependence of the density contrast k (equation J3|) on 
radius also has an influence on the details of breakout. When 
the dissociation energy e is held fixed, k increases toward 
larger radius. This has the effect of creating a dense layer 
of fluid below the shock when the shock has reached a radius 
where e ~ vj f /2. In spherical symmetry, the breakout of the 
shock is then impeded by this layer, which cannot exchange 
position with the lighter material below it. It can happen that 
the energy in the expanding region is no longer able to sustain 
the heavier material above, and the shock collapses, as shown 
in Figure|4]for the constant-e run with H = 1.08// cr . This ob- 
struction is avoided when statistical equilibrium between n, 
p, and a is maintained below the shock, because e/v\ and 
K both decrease gradually as the shock expands to distances 
much larger than r s o (Figure[TJ. This limit to the shock expan- 
sion does not occur in two dimensions, as the superposition of 
dense fluid over lighter fluid is Rayleigh-Taylor unstable on 
the dynamical time tfjx>. 

4. TWO-DIMENSIONAL SIMULATIONS 

Extending the flow calculation to two dimensions reveals 
some subtle patterns of behavior. The time evolution of the 
shock is shown in the left panel of Figure|6]for a range of heat- 
ing rates near the threshold for explosion. In contrast with the 
one-dimensional runs, the breakout of the shock looks similar 
in models with constant dissociation energy and with NSE be- 
tween n, p, and a below the shock. Both types of models are 



subject to buoyancy-driven instabilities, which allow cold ma- 
terial below the shock to interchange position with hotter ma- 
terial within the gain region. As a result, the shock is highly 
asymmetric at breakout in both cases. In {(5]we compare the 
critical heating rate for explosion in one- and two-dimensional 
runs, and examine how it is influenced by a-particle recombi- 
nation. 

Around the threshold for explosion, all of our runs develop 
vigorous convective motions before the SASI has a chance to 
undergo even a few oscillations. At high heating rates, we find 
that the convective instability is driven by the negative entropy 
gradient within the layer of maximal neutrino heating. In 
non-exploding runs, the shock settles to a quasi-equilibrium 
state with oscillations taking place over a ran ge of angu- 
lar (Legendre) index £, as previ ously seen b y Ohnish iet alj 
(2006). IScheck et all (120081) . and lMurphy & Burrowsl (12008). 
The amplitude of the I = 1 and 2 modes remains small until 
the heating parameter H has begun to exceed about one half 
the critical value for an explosion. The competition between 
SASI growth and convective instability is examined in detail 
in m 

One gains considerable insight into the mechanism driving 
shock breakout by examining the distribution of Bernoulli pa- 
rameter (equation G)) in the shocked fluid. We first consider 
the NSE runs with r s o = 50 km and 125 km, with the heat- 
ing parameter H just above the threshold for an explosion. 
Two snapshots from each of these runs are shown in Figure|7] 
In the first case, the initial equilibrium shock radius is only 
~ r Q /4 km, and a-particles are essentially absent below the 
shock. In the second, the shock starts at ~ 2r a /3 and X a ~ 0.5 
initially in the postshock flow. 

Large deformations of the outer shock are caused by con- 
vective plumes that carry positive energy. Strong neutrino 
heating is generally concentrated within an inner zone where 
the material is gravitationally bound (b < 0). The degree of 
symmetry of this bound material depends on the a-particle 
abundance. In the r s o =125 km run, it is spherically strati- 
fied and the material with b > is generally excluded from it. 
Strong recombination heating is present both below and above 
the surface where b ~ 0, indicating that it is mainly responsi- 
ble for imparting positive energy to the shocked material. The 
mean shock radius expands by a factor ^2.5 between the up- 
per two frames in Figure |7l but the growth in the volume of 
positive-energy material is not accompanied by a significant 
expansion of the inner bound region, whose outer radius re- 
mains fixed at r ~ r a . 

This segregation of bound from unbound material is broken 
when the shock is more compact initially, as is seen in the 
lower two panels of Figure [7] A single dominant accretion 
plume is continuously present, which funnels cold and dense 
material into the zone of strong neutrino heating. Alpha- 
particles are present only well outside the boundary between 
b < and b > 0. The competition between a -part icle and 
neutrino heating is discussed in more detail in 34.11 and the 
influence of a-particles on the threshold heating rate for an 
explosion is examined in $5] 

The accumulation of a bubble of hot material right behind 
the shock is a consequence of the balance of the buoyancy 
force acting within the bubble, and the ram pres sure of the 
presh ock material. The ratio of force densities is dTho mnson 
120001) 

^buoy / P-Pbubble ^ ( 2GM/r s \ 

-5 — i AiZbubble, (15) 
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FIG. 7. — Snapshots of two separate models, with heating parameter //just above the threshold for an explosion, and initial shock radius either well inside r a 
(r s o = 50 km, without heating) or close to r a (r s o = 125 km, without heating). Within each panel, the top figure displays Bernoulli parameter b; the middle figure 
the rate of change of nuclear energy generation; and the bottom figure the net rate of neutrino heating. Top left: the early development of an asymmetric plume 
with positive b; top right: the same run just before the shock hits the outer boundary. In this r,o = 125 km run, the heating by a-particle recombination is enhanced 
with respect to neutrino heating due to the large X a in the initial stationary model. The central zone with b < maintains a nearly spherical boundary near the 
radius r a (~ 2.0f 5 o), and recombination heating straddles this boundary. Bottom left: a-particles begin to form as the shock approaches r a in the r,o = 50 km 
run, but neutrino heating remains much stronger than recombination heating. Bottom right: the same run just before the shock hits th e out er boundary. When the 
shock starts off well inside r a , neutrino heating dominates the initial expansion, and material with b > forms well inside r a (see 34.3) . Animations showing 
the evolution of these two configurations are available in the online version of the article. 
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FIG. 8. — Normalized pressure gradient (r//>)|Vp| showing the secondary 
shock structure during breakout. The top model is e = 0. 15vj| „ and the bottom 
NSE with r s o = 75 km. Both have heating rates just above the threshold for 
an explosion. An animation showing the time evolution is available in the 
online version of the article. 
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FIG. 9. — Top panels: rate of release of specific nuclear binding energy 
de nuc /d;. Bottom panels: mass fraction of a-particles X a . We show two 
instants in the exploding NSE run with r s o = 75 km and H = 1.02// cr . The 
shock contour is approximated by the white line which marks Xq = 90%. 



where r s is the shock radius, v r is the ambient radial flow 
speed, p is the ambient density, ^bubble the density of the bub- 
ble, and A^bubbie is its angular size. A low-density bubble 
(/O-pbubbie ~ p) can resist being entrained by the convective 
flow once it grows to a size Afibubbie ~ -Ml oa Sr, where Ai CO n 
is the convective Mach number. On the other hand, the bubble 
must attain a much larger angular size ASlbubbie ~ 1 Sr if the 
buoyancy force is to overcome the upstream ram (|v r | ~ Vff) 
and force a significant expansion of the shock surface. Fig- 
ure [7] shows that the extent of the shock expansion is indeed 
correlated with the angular width of the region where hot ma- 
terial accumulates. 

Another interesting feature of Figure [7] is the presence 
of secondary shocks, which are triggered once the outer 
shock becomes significantly non-spherical. Their locations 
are marked by discrete jumps in the rate of recombination 
heating. Secondary shocks are also prevalent throughout the 
nonlinear phase in the constant-e models. Figure [8] shows the 
normalized pressure gradient (r/p)\Vp\ for collapse models 
of both types, when H is just above threshold for an explo- 
sion (right before the shock hits the outer boundary of the 
simulation volume). The online version of the article contains 
an animated version of Figure [8] showing the complete evolu- 
tion. In both cases, secondary shocks extend over the whole 
postshock domain, signaling the dissipation of supersonic tur- 
bulence which is stirred by accretion plumes that penetrate 
into the gain region. 

4.1. Distribution of Alpha Particle Recombination Heating 

Heat input by neutrino absorption and by a-particle recom- 
bination have very different distributions within the shocked 
fluid: strong neutrino heating is concentrated inside r s o, 
whereas recombination heating of a comparable amplitude 
is distributed throughout the settling flow. Strong recombi- 
nation heating quite naturally extends below the zone where 
a-particles are present in significant numbers, as is seen in 
Figure [9] The first and third panels of this figure depict the 
pre-explosion steady state of the r s o = 75 km model with 
H = 1 .02H CI , while the second and fourth panels show the 
last time before the shock hits the outer boundary. At the 
latter time, one sees that the strongest recombination heating 
is concentrated in a layer where X a < 0.5, at the base of the 
extended a-rich plumes. Just as in the one-dimensional sim- 
ulations (e.g. Figure 0, X a approaches unity during shock 
breakout. 

The relative strength of neutrino heating and recombina- 
tion heating depends on the initial radius of the shock, and 
on the Bernoulli parameter of the postshock material. Fig- 
ure [10] separates out cooling by a-particle dissociation from 
heating by recombination and neutrino radiation during the 
pre-explosion quasi-steady state (leftmost panels), at the onset 
of explosion (second panel left to right), and during breakout 
(two rightmost panels). See Figure [6] for comparison. The 
colored curves show the positive, negative, and net contri- 
butions from nuclear energy generation. The sharp negative 
spike near b = represents a-particle dissociation in fresh, 
cold downflows. The formation of material with b > is pri- 
marily due to a-particle recombination in the r s o =125 km 
run. As the initial radius of the shock is reduced with respect 
to r a , neutrino heating makes a proportionately larger contri- 
bution near breakout. 

The strength of the boost given to the shock by recombina- 
tion heating can be gauged by comparing de nuc /df to the adi- 
abatic rate of change w a d of the enthalpy of the flow (equation 
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FIG. 10. — Heating rate of material, as distributed with respect to Bernoulli parameter b. This illustrates the relative importance of neutrino heating and 
nuclear dissociation/recombination in hot and cold parts of the flow. We restrict attention to material in the gain region (defined by J£n > _S?c) in me three 
two-dimensional NSE runs just above the threshold for explosion. Four snapshots are shown: the pre-explosion quasi-steady state (leftmost), onset of explosion 
(second from left to right), and breakout (third and fourth). Black curves: net heating rate resulting from neutrino absorption and emission. Red/green curves: 
heating/cooling rate by a-particle recombination and dissociation in material with de nuc /d/ > and de nuc /dt < 0, respectively. Blue curves: net heating/cooling 
rate due to changing a-particle abundance. The sharp negative spike near b = represents a-particle dissociation in fresh, cold downflows. The formation of 
material with b > is primarily due to a-particle recombination in the r s g = 125 km run. As the initial radius of the shock is reduced with respect to r a , neutrino 
heating makes a proportionately larger contribution near breakout. 



Ifl4l ). Figure [TTI shows the result for all three NSE sequences 
with //just above H a . In all cases, de nuc /df ~ w a d in various 
parts of the shocked fluid once the shock extends beyond a 
radius ~ r a . Most of the heat input by recombination is con- 
centrated where X a ~ — 0.5, just as in spherical symmetry. 

Histograms of de nuc /df versus b and X a are shown in the 
right panel of Figure [12] The rapid dissociation of a-particles 
in fresh downflows is represented by the long tail toward large 
negative values of de nuc /df, showing that the overall contribu- 
tion of nuclear energy generation is negative. The a-particle 
concentration is very stratified, with higher X a occurring at 
larger radius. Most of the mass with positive Bernoulli pa- 
rameter is located at large radii. 



4.2. Residency Time 

A long residency time of material in the gain region is com- 
monly viewed as a key ingredient in a successful neutrino- 
driven expl osion. To calculate f res , we use the method de- 
scribed in 32.1.31 we first assign a unique "fluid time" tf 
(equation |[T2"1 ) to each infalling radial mass shell in the sim- 
ulation, which is effectively the time at which it passes the 



shock. We then define 6 the residency time of the fluid as 

t KS = t-t F , (16) 

where t is the p resent time. A related meth od (tracer parti- 
cles) is used by iMurphv & Burrowsl ([2008 ) to calculate the 
residency time in collapse simulations with a more realistic 
EOS. 

As material with positive Bernoulli parameter accumulates 
below the shock, we indeed find that its f res grows larger. 
The shock starts running outward if the energy of this un- 
bound material grows on a timescale shorter than the con- 
vective time. The right panel of Figure [6] shows the charac- 
teristic expansion time of the shock f exp ~ r s max /|dr s inax /df | 
(as measured at its outermost radius), alongside (? res )voi (as 
measured within the material comprising the upper part of the 
residency time distribution). The final breakout of the shock 
seen in the left panel of Figure[6]corresponds to the time when 
fexp ~ (fies)voi- This lengthening of the mean residency time 
can largely be ascribed to the increased dynamical time of 
the expanding shock. What changes most dramatically during 
breakout is the ratio of the expansion time to the dynamical 

6 Since tf is defined in terms of the angle-averaged radius of the shock, 
there is a modest error in tf due to non-radial deformations of the shock. 
Given the lack of substantial large-scale mixing between the single accre- 
tion funnel and convective cells, this prescription serves well as a tracer of 
different fluid populations. 
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FIG. 11. — Ratio of de nl | C /df (rate of release of nuclear binding ener: 
equation (7J) to w a( ] (adiabatic rate of change of the enthalpy, equation 
NSE models shown have H just above H CI . 

time. 

The breakout is a bit more gradual in the r s0 =125 km 
model with heating just above threshold for an explosion 
(H = \.QAH a ; see animated version of Figure |7^,b in the on- 
line material). In this case, the expansion time of the shock 
remains somewhat longer than the residency time of material 
below the shock, which implies that the breakout depends on 
the continuing release of nuclear binding energy. 

Note that large changes in the distribution of f res are concen- 
trated in regions of positive b. In the left panel of FigurjT2l 
we plot the distribution of b and f res in the e = 0.15v^ run that 
is just above the threshold for an explosion. Regions with 
small or negative residency time represent freshly injected 
fluid. The distribution is stratified in b around f res ~ 40fffo 
(which is app roximately an overturn period of a convective 
cell, see 34. 3 1 and Figure [13}. Material with more negative b 
resides on average at a smaller radius. Fluid with a longer res- 
idency time has mostly positive b, corresponding to material 
transported upwards by convective cells. 

It is also apparent from the right panel of Figure[T2]that ma- 
terial with a longer residency time tends to have lower X a , as 
is expected because it also tends to have a higher temperature. 

4.3. Heat Engine in a Two -Dimensional Explosion 

Fresh material that is accreted through an oblique shock 
has a relatively low entropy, but once it reaches the base of 
the gain region it is exposed to an intense flux of electron- 
type neutrinos. Some of this heated material rises buoyantly, 
and forces an overturn of the fluid below the shock. Material 
with a longer residency time may t herefore undergo multiple 
episodes of heating. On this basis. iHerant et alj (fl 99211 19941) 
suggested that a convective flow would mediate a heat engine 



below the shock that would drive a secular increase in the en- 
ergy of the shocked fluid. 

We now investigate whether a heat engine operates in our 
simulations, and how it depends on the heating parameter 
H. We focus on a model with a constant nuclear dissocia- 
tion energy, e = 0.15vj| . In this class of models, the infalling 
heavy nuclei are completely broken up below the shock, and 
no heating by the reassembly of a-particles is allowed. As 
a first step, we average the convective flow over windows of 
width 50fffo, which de-emphasizes short term fluctuations in 
the averaged velocity field (v),. Figure[T3lshows (v), and (b), 
(equation [0) at four different times in the run with H just 
above the threshold for explosion (H = 1.02// cr ). The radius 
of maximum heating (r ~ 0.66r s o) coincides with the lower 
boundary of the convective cells, across which material flows 
horizontally. The overturn period in these large scale cells is 
~ 40—50%), as found by integrating streamlines of the mean 
flow (an example is shown on the lower-left panel of Fig- 
ure [T3l . Heated fluid accumulates in the region in between 
the top of convective cells and the shock. A strong defor- 
mation of the shock allows a plume of fresh material to de- 
scend diagonally between the convective cells. The tilt of this 
cold downflow intermittently flips in sign, and the averaged 
circulation pattern typically has an "oo" shape. The heating 
of fluid parcels in the two hemispheres is also intermittent, 
and sometimes two circulation flows are established simulta- 
neously, thereby causing a bipolar expansion of the shock. 

We have identified a useful figure of merit which connects 
a secular increase in the shock radius to the strength of neu- 
trino heating at the base of the gain region. Figure [14] shows 
the absolute value of the Bernoulli parameter b at the radius 
of maximum heating (re.max — 0.66r s o) as a function of po- 
lar angle 6. In the top panel, the four sets of thin solid lines 
correspond to the four snapshots of FigureQj] and the bottom 
panel shows the analogous results for a non-exploding run. 
Overplotted as thick solid lines is the quantity 



e = 



>'H.i 



\{ve)t,8* 



(17) 



This measures the specific energy that is absorbed from neu- 
trinos by the material that flows laterally along the lower 
boundary of the convective cells. In equation dl7} . the an- 
gular average of the heating rate and meridional velocity is 
restricted to a single convective cell. 7 

In non-exploding models, the circulation in the gain region 
settles to a quasi-steady state, with no net amplification of the 
mass in material with positive b. The heat absorbed at the base 
of the convective cells is of the same order of the Bernoulli 
parameter of the fluid, that is, < \b\ . In an exploding model, 
9 will often exceed \b\ by a factor 2-3. As is shown in the 
upper panel of Figure [14] grows with time as the system 
approaches the explosion. 

The stability of the averaged flow pattern appears to be, in 
part, an artifact of the axisymmetry of the flow. This im- 
poses strong restrictions on the motion of convective cells, 
causing vorticity to accumulate on the largest spatial scales. 

7 To identify the range of angles comprising the lower boundary of a con- 
vective cell (r = rn.max), we first average (vg)t over all polar angles, and then 
define a single cell as a zone where |(v,-) f | < | (ve) f ,e | and vg maintains a 
constant sign. Once the convective cells have been so identified, the angu- 
lar avera ge is repeated within each cell. The quantity \{vg), g t | appearing in 
equation (TTJ represents this more restrict ed a verage, which typically covers 
~ 1 rad in the polar direction (e.g., Figure [T3l . 
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FIG. 12. — Left panel: Histogram of Bernoulli parameter b and residency time f res in the exploding run with e = 0.15i»f f0 (see also Figure [T3V The colors label 
the mass-weighted radius, and we include all material experiencing a net excess of neutrino heating over cooling. Right panel: Histogram of Bernoulli parameter 
b and residency time ( lcs versus a-particle mass fraction X a and rate of release of specific nuclear binding energy de nuc /df , in the exploding run with r s a = 75 km 
(Figure [9). 
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FIG. 13. — Four snapshots of the exploding model with e = 0.15Vf f0 and H = 1.02// cr - The Bernoulli parameter (color map) and the velocity field (white 
arrows) are averaged over intervals of duration 50%o. The thick white contours show the surface with 50% mass fraction in heavy nuclei (time-averaged). The 
yellow curve in the lower-left panel shows the result of integrating a streamline of this time-averaged velocity field, starting from a point just above the radius of 
maximum heating. The curve performs an overturn after ~ 45%0i an d takes an extra ~ 8%q to reach the inner boundary. 
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FIG. 14. — Bernoulli parameter (thin solid line s) and specific energy ab- 
sorbed during lateral advection (equation 1171 . thick lines) at the radius 
of maximum heating r^max- These quantities are averaged over intervals of 
duration 50fffo- The upper panel shows the exploding run with e = 0.15v'f f0 
and H = l.Q2Hcr (the same as in Figure \l3\ . and the lower panel shows the 
non-exploding run with e = 0.15Vf f0 and H = Q.9Lffcr. Dotted lines show the 
angular boundaries of convective cells. 

Our observation that the bulk of the neutrino heating takes 
place within horizontal flows suggests that the ratio of heat- 
ing timescale to radial advection time in the gain layer may 
be a less precise diagnostic of the conditions for explosion 
in two dimensions: the horizontal convective velocity is typ- 
ically low compared with the downward velocity of the main 
accretion plume. We do observe that the main accretion plume 
becomes strongly distorted near the threshold for an explo- 
sion, so that a significant fraction of the plume material en- 
ters one of the convection cells. This effectively decreases the 
amount of material that accretes to the protoneutron star and 
thus increases the overall advection timescale across the gain 
region. 

5. CRITICAL HEATING RATE FOR EXPLOSION 

An explosion occurs when the heating parameter// is raised 
above a critical value 8 H a . We now explore how H cl de- 
pends on the details of the EOS and the initial radius of 
the shock. One can express H simply in terms of the ra- 
tio of the heating rate (47rr 3 Jz?//) to the accretion luminosity 
(GMM/r), in the idealized (but unrealistic) case where the 
flow is composed only of free nucleons and moves hyperson- 
ically. Then this ratio depends on H but not on the accretion 
rate M = 47rr 2 p(r)|vff(r)|. The precise value of the reference 

8 Our method for determining H a is discussed in j|3.1l 
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FIG. 15. — Critical heating parameter H a that yields an explosion, for all 
the model sequences explored in this paper (TableQ}. The abscissa is the ratio 
of e to v\ in the initial flow configuration (vi being the radial flow velocity just 
upstream of the shock). Error bars show the separation between exploding 
and non-exploding models, with the points marking the average. 

radius is unimportant; we choose r,o, the shock radius in the 
time-independent, spherical flow solution at H = 0. Then 



47rr s0 3 JfH[/3i(^o)] 



GMM / r S Q 



X a =0; M=oo 



2H 



(18) 



This quantity is ~ 10~ 3 - 10~ 2 in the models we examine, 
which are below or near the threshold for explosion. 

Note that the cooling in our model is concentrated at the 
base of the settling flow. As a result, the width of the gain 
region (relative to the shock radius) does not change signifi- 
cantly between different models. The critical heating param- 
eter is therefore only indirectly related to the amplitude of 
the cooling function through the structure of the settling flow 
below the shock. Our purpose here is to explore how the crit- 
ical heating rate depends on the strength of the gravitational 
binding of the shocked fluid to the collapsed core, and on the 
abundance of alpha particles. 

Figure Q3] displays H a for all of our model sequences. The 
abscissa is e/v 2 , where e is the nuclear dissociation energy 
and Vi is the flow speed upstream of the shock in the initial 
configuration (that is, in the time-independent, spherical flow 
solution). In the case of the NSE equation of state, this quan- 
tity can be translated into an initial value of the shock radius 
using Figure [TJ (Note that e/v\ has a weak dependence on r s 
in the NSE sequence.) 

A few interesting features of Figure [15] deserve comment. 
First, a comparison with TableQ]shows that the critical heating 
rate for explosion is ~ 50-70% of the maximum heating rate 
for which a steady-state flow solution can be found. The max- 
imal heating parameter Z/ Ste ady for a stead y flow corresponds 
directl y to the one first determined by lBurrows"& Goshv 
(1993)" using a more realistic EOS. Note also that the values of 
H cr in the one- and two-dimensional models are much closer 
to each other than they are to //steady This result is perhaps 
not surprising, given that the explosion is not immediate, but 
is approached through a series of transient fluid motions. 

Second, H a is lower when NSE between n, p and a is main- 
tained below the shock. In this case, the dissociation energy at 
the shock is not fixed, but is (roughly) inversely proportional 
to radius. However, the difference between the NSE mod- 
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FIG. 16. — Critical heating parameter H a that yields an explosion, for the 
runs that include a-particles in the EOS. The abscissa is the ratio of t he in itial 
shock radius r s to r a . Error bars have the same meaning as in Figure [T31 The 
critical heating parameter (a close analog of L v ) decreases substantially with 
increasing shock radius. The differences in Ha between the one- and two- 
dimensional models also decreases. 

els and the constant-e models is only ~ 10% in H a when the 
shock starts out well below r a (equation HI). An explosion is 
significantly easier when the fluid below the shock starts out 
with a significant population of a particles, as in the models 
with r s o = 125 km. 

Third, H a tends to decrease with increasing e/v\: a slightly 
lower heating rate per unit mass is required to explode a 
flow with a larger density contrast k across the shock. Be- 
cause almost all the gravitating mass is in the collapsed core, 
the gravitational binding energy of the gain region is ap- 
proximately proportional to k, whereas the net heat absorbed 
over the advection time is a stronger function of density, 
fadv / (=% - -2c)d 3 '" oc k 2 . (One factor of k comes from the 
advection time f a d v as given by equation 11131 . and the other 
from the density dependence of J?h-) For example, Table [TJ 
shows that k is ~ 1.6 times larger for e/vg = 0.2 than for 
e /vgQ = 0.1, and that H cl is smaller by the inverse of the same 
factor. 

Fourth, the two-dimensional runs all require less heating 
than their spherically symmetric counterparts to explode. A 
major reason for this is that all two-dimensional configura- 
tions explode along one or both poles (see Figs. [TT1 and [TJll. 
so that less material must be lifted through the gravitational 
field than in a fully spherical explosion. We have found that 
the precise value of the difference between the critical heating 
rate in the one- and two-dimensional explosions depends on 
the choice of r» / r s o, and therefore on the normalization C of 
the cooli ng function. The fact that w e find a smaller differ- 
ence than Mur phy & Burrowsl (120081) may be a consequence 
of our simpler cooling function and equation of state. 

The critical heating rate depends in an interesting way on 
the starting radius of the shock, in a way that points to the 
recombination of a-particles as an important last step in the 
transition to an explosion. Figure [16] shows that H a in the 
NSE models grows rapidly as the initial shock radius 9 r s is 
pushed inside r a . Here we normalize the heating parameter at 
a fixed physical radius, namely r a . Translated into the context 

9 Note that r, is the shock radius in the time-independent flow solution. 
For a fixed cooling function, r, is a monotonically increasing function of H, 
and equals r s Q at H = 0. 
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FIG. 17. — The development of a convective instability is strongly limited 
when the parameter \ < 3. These panels show snapshots of entropy (normal- 
ized to initial postshock value) in a NSE run (r 5 o = 75 km) with two different 
heating rates. When \ = 2.1, convective cells of a limited extent are trig- 
gered in the layer where the net heating rate is strongest, but they do not 
propagate into the upper parts of the gain region. Convection becomes much 
more vigorous and widespread when x = 5.5. Note that both of these models 
are non-exploding. An animation showing the time evolution of these two 
configurations is available in the online version of the article. 

of a realistic core collapse, this means that the critical neu- 
trino luminosity for an explosion decreases with increasing 
shock radius. The radius of the stalled shock depends, in turn, 
on the EOS above nuclear matter density: iMarek & Jankal 
(2009) find that a softer EOS corresponds to a larger shock 
radius, mainly due to the higher accretion luminosity onto the 
neutronized core. Here we have subsumed this uncertainty in 
the high-density EOS into a single free parameter, the ratio 
r s o/r a . Hydrodynamic instabilities are effective at driving an 
explosion to the extent that they push the shock radius close to 
r a ; beyond this point, the remainder of the work on the flow 
is done largely by a-particle recombination. 

One also notices from Figure[T6]that the difference between 
H CI in one and two dimensions depends on the starting radius 
of the shock. The closer r s o is to r a , the weaker the depen- 
dence of the critical heating rate on the dimensionality of the 
flow. 

6. CONVECTION AND THE SASI 

Overturns of the fluid below the shock can be triggered in 
two distinct ways: through the development of Ledoux con- 
vection in the presence of a strong negative entropy gradient, 
or via the non-linear development of the SASI (a linear feed- 
back between ingoing entropy and vortex waves, and an out- 
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FIG . 18. — Amplitude (r.m.s.) of £ = 0, 1 , 2 modes of the shock in two model 
sequences with varying heating parameter H. Stars indicate exploding runs. 
We show the r.m.s. fluctuation of the difference between the instantaneous 
Legendre coefficient ag and a running average (ae)sot that is computed over 
a window of width 50%o (see text). This subtracts the secular movement of 
the shock in runs that are close to or above the threshold for explosion. Note 
that the amplitude is measured in absolute units (r s o). 

going sound wave). We now show that the amplitude of the 
dipolar mode that is excited in the shock is strongly tied to 
the level of neutrino heating, and so thermal forcing plays a 
crucial role in maintaining the oscillation. To a certain ex- 
tent, this distinction is of secondary importance, in the sense 
that memory of the linear phase of the instability is lost once 
the inflow of fresh material below the shock bifurcates from 
older shocked fluid. Nonetheless, the origin of the convective 
motions does have implications for the stability of £ = 1 and 
2 modes in three-dimensional simulations: one expects that 
large scale oscillations will change shape and direction more 
rapidly if they are triggered primarily by neutrino heating. 

We can ask whether a heating parameter H that yields an ex- 
plosion will also form an unstable entropy gradient below the 
shock. Convection develops through a competition between 
in ward advection and n eutrino heating. A detailed analysis 
by Foglizz o et alj (12006) shows that the parameter 
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(19) 



must exceed a critical value ~ 3 for unstable convective 
plumes to grow before being advected downward through the 
gain region below the shock. Here wbv is the Brunt-Vaisala 
frequency. 

Using our initial flow models, we can translate H into a 
value for \, and find (Table [TJ that typically \ ~ 5-20 at the 
threshold for a neutrino-driven explosion. The implication for 
convection below the shock is illustrated in Figure [17] which 
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FIG. 19. — Principal shock Legendre coefficient and the radial median of 
the angle-averaged entropy gradient for the e = 0.15vf f0 model and two dif- 
ferent heating rates, corresponding to \ = 1 .5 (H = 0.002 r s o Vj™, blue curves) 
and x = 3.9 (H = 0.004r j0 v| , red curves). Both runs are below the threshold 
for an explosion, but vigorous convection is established throughout the gain 
region in the run with the higher heating rate. Top (bottom) two panels: seed 
perturbation is a shell with I = 1 (£ = 2) density profile. A running average of 
the entropy gradient (temporal width 20fffo) appears as thick solid lines. 

shows two snapshots for models with \ = 2.\ and 5.5, nei- 
ther of which explodes. At the lower heating rate, the time 
required for convection to develop depends on the strength of 
the seed perturbation, whereas at the higher heating rate con- 
vective overturns develop rapidly within the layer of strong 
neutrino heating and spread throughout the post-shock region 
over a few dozen dynamical times. (The figure shows the re- 
sult in the case where the seed perturbation is dominated by a 
small spherical startup error in the initial model.) 

We conclude that, near the threshold for a neutrino-driven 
explosion and for our given set of physical assumptions, 
convection is driven primarily through the development of 
a strong, negative entropy gradient within the gain region, 
rather than through the non-linear development of SASI 
modes. The growth of the SASI requires at least a few oscil- 
lations, each with a period comparable to the advection time. 
The SASI is therefore subdominan t when y > 3. It is worth 
comparing this with the results of Sch eck et al.1 (120081) . who 
forced the inner boundary of the simulation volume to move 
inward to model core contraction, thereby generating large 
advection velocities. The net result was that the flow barely 
reached x — 3 in exploding configurations. While this effect 
may be important for relatively prompt explosions, it should 
be kept in mind that the rate of contraction of the neutri- 
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FIG. 20. — Same as Figure [191 but for NSE models with a-particles and 
r s0 = 75 km. 

nosphere has slow ed subtantially a few h undred milliseconds 
after core bounce. iMarek & Janka (2009) found evidence for 
the delayed explosion of a 15M Q progenitor around 600 ms 
after core bounce, for which the effect of core contraction is 
not likely to dominate the dynamics. 

In Paper I we considered the non-linear, saturated state of 
the SASI in the absence of neutrino heating, and showed that 
the amplitude of the shock oscillations drops significantly as 
the dissociation energy e is increased. We now explore how 
the r.m.s. amplitude of the shock oscillations correlates with 
the strength of heating. To eliminate the effect of secular 
shock motions around or above the threshold for explosion, 
we first calculate the running average {ai)sot of the shockLeg- 
endre coefficients ai over an interval 50%o, and then calculate 
the r.m.s. of at = at - (ai)sot over the duration of each simu- 
lation. The result is plotted in Figure[l8]as a function of H for 
two model sequences (e = 0.15Vg ar >d NSE with r s o = 75 km). 

There is a clear trend of increasing £=1,2 mode ampli- 
tude with increasing heating rate. This confirms our previ- 
ous suggestion that large-amplitude shock oscillations require 
strong heating when the dissociation energy behind the shock 
exceeds ~ 0.15vg . The models with H = reveal a slight ex- 
ception to the overall trend: the r.m.s. amplitude of the shock 
oscillations appears larger than it does in models with small 
but finite heating rate, because the oscillations are strongly in- 
termittent at H = (see Paper I). The shock oscillations grow 
much stronger just below the threshold for explosion (explod- 
ing runs are marked by stars), above which they seem to satu- 
rate. Note also that their amplitudes do not vary much with the 
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FIG. 21. — Histogram of vorticity vs. Mach number in the gain region, 
weighted by mass. Shown are three different instants in the evolution of the 
run with e = 0. 15vg , £ = 1 see d pe rturbation, and \ = 3.9 (corresponding to 
the upper two panels in Figure [l9l . The distribution broadens once convec- 
tion is fully developed, but before the dipolar shock mode shows significant 
growth. The vertical lines show the vorticity of a convective flow with period 
equal to (solid) the mean radial advection time and (dashed) the period of a 
lateral sound wave at the midpoint between r„ and r s . An animated version 
of this figure is available in the online version of the article. 

choice of dissociation model. The r.m.s. amplitudes relative 
to the running average of ao at the threshold for explosion are 
{5%, 12%, 8%} for the 1= 1,2,3 modes in the e/v| = 0.15 
sequence, and {6%, 12%, 7%} in the NSE r s o = 75 km se- 
quence. 

We have performed an additional sequence of runs in which 
we drop an overdense shell with a given Legendre index I 
through the shock (see also Paper I). This has the effect of 
selectively triggering individual SASI modes. Figures [19] and 
l20l display the Legendre coefficients of the shock alongside 
the angle-averaged entropy gradient. We find that the ampli- 
tude of the £ = 1 and 2 modes is strongly tied to the strength 
of convection. For both types of dissociation models, convec- 
tion is quenched by the accretion flow when \ < 3: it grows 
intermittently in strength, but never reaches large enough am- 
plitudes to significantly distort the shock surface. As a con- 
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sequence, the entropy gradient remains shallow and negative 
most of the time. Coherent shock oscillations are also seen, 
but they have a low amplitude due to the large dissociation 
energy. 

Convection grows much more rapidly when \ > 3, dis- 
torting the shock surface before the SASI has the chance to 
execute a few oscillations. In this case, the entropy gradi- 
ent is initially more negative, but quickly flattens. Indeed, 
the £ = 1,2 amplitudes only become large in models where 
neutrino-driven convection is strong enough to flatten the en- 
tropy gradient. 

Another way of seeing that convection is the forcing agent 
behind shock oscillations when x > 3 is to analyze the dis- 
tribution of vorticity in the gain region. Figure [21] shows a 
histogram of vorticity vs. Mach number at three different in- 
stants in the evolution of the run with e = 0.15vj f0 , £ = 1 seed 
perturbation, and \ = 3.9 (upper two panels in Figure [T91. At 
t = 20fffo, convection is just getting started and the vortical 
motions are restricted to Mach numbers < 0.3. However, by 
t = 35fffo the Mach number distribution extends up to M. > 0.5 
and has almost reached its asymptotic form (t = 60%)), at 
the same time that convection has filled the region below the 
shock. The dipolar mode of the shock develops a large ampli- 
tude only after this fully developed convective state has been 
reached. The convective rolls are a source of acoustic radia- 
tion (e.g. lGoldreich & Kumar 1988), which will drive a dipo- 
lar oscillation of the shock if the overturn frequency is compa- 
rable to the frequency of the £=\ mode, | V x v ~2x lir/t^. 
This zone is marked by the vertical solid lines in Figure I2T1 
and indeed encompasses most of the mass. 

7. SUMMARY 

We have investigated the effects of a-particle recombi- 
nation and neutrino heating on the hydrodynamics of core- 
collapse supernovae, when the heating rate is pushed high 
enough to reach the threshold for an explosion. The effect of 
dimensionality has been probed by comparing one- and two- 
dimensional time-dependent hydrodynamic calculations. Our 
main results can be summarized as follows: 

1 . - The critical heating parameter that yields an explosion 
depends sensitively on the starting position of the shock rel- 
ative to r a . This means that the critical neutrino luminosity 
depends sensitively on the stall radius of the shock and, in 
turn, on the core structure of the progenitor star and the den- 
sity profile in the forming neutron star. Within the framework 
explored in this paper, we find two extreme types of explo- 
sion. In the first, neutrino heating does most of the work, 
with a significant final boost from a-particle recombination. 
In the second, neutrino heating is generally less important at 
promoting material below the shock to positive energies. 

2. - During the final stages of an explosion, the heat released 
by a-particle recombination is comparable to the work done 
by adiabatic expansion. This heat is concentrated in material 
that has previously been heated by neutrinos. Significantly 
more energy is lost through a-particle dissociation in fresh 
downflows, so that nuclear dissociation remains on balance 
an energy sink within the accretion flow. 

3. - The large-amplitude oscillations that are seen in one- 
dimensional runs near an explosion are the consequence of 
the £ = SASI as modified by heating. In contrast with the 
1=1,2 modes of a laminar accretion flow, the period of these 
oscillations is close to twice the post-shock advection time. 
The critical heating rate for an explosion (assuming constant 



mass accretion rate, neutrino luminosity, and inner boundary) 
corresponds to neutral stability for the £ = mode. 

4. - The critical heating parameter H a for an explosion is 
generally lower in two dimensions than in one, but the dif- 
ference becomes smaller as the starting radius of the shock 
approaches r a . The difference depends somewhat on the ratio 
r*/r s o and thus on the cooling efficiency and equation of state. 

5. - Non-spherical deformations of the shock are tied to the 
formation of large-scale plumes of material with positive en- 
ergy. Our two-dimensional explosions with a super-critical 
heating rate involve a large-scale convective instability that 
relies on the accumulation of vorticity on the largest spa- 
tial scales. Volume-filling convective cells are apparent in a 
time-averaged sense. Transient heating events create positive- 
energy material that accumulates in between the convective 
cells and the shock. A significant fraction of the heating oc- 
curs in horizontal flows at the base of the convective cells, 
which are fed by a dominant equatorial accretion plume. If 
the heating parameter is large enough, this results in an am- 
plifying cycle and explosion. 

6. - The amplitude of the £ = 1 and 2 modes correlates 
strongly with the value of the heating parameter, and is 
coupled to the appearance of vigorous neutrino-driven con- 
vection below the sh ock. In agreement wit h the work of 
iFoglizzo et all d2006l) and IScheck et al.1 d2008l) . we find that 
X~3 marks the transition from a strong linear instability in a 
nearly laminar flow below the shock, to a volume-filling con- 
vective instability. In all of our simulations, the threshold for 
explosion lies well within the latter regime. This highlights 
a basic difference between one- and two-dimensional explo- 
sions: the mechanism is fundamentally non-linear in two di- 
mensions. 

7. - We have explored essentially one ratio of cooling radius 
to shock radius, namely r*/r s o = 0.4 at zero heating (corre- 
sponding to r*/r s ~ 0.2 near the threshold for an explosion). 
The growth of the £ = 1 SASI mode is strongest for this partic- 
ular aspect ratio when e = (see Figure 12 of Paper I). As dis- 
sociation is introduced into the flow, we found that the peak 
growth rate moves to larger values of r*/r s o. On the other 
hand, detailed collapse calculations indicate ratios of neutri- 
nosphere radius to shock radius that are ev en smaller than 
~ 0.2 following ~ 100 ms after collapse (e.g. lMarek & Jankal 
120091) . We conclude that the I = 1 SASI mode is not being 
artificially suppressed by our choice of initial shock size. 

8. - Vortical motions with a Mach number ~ 0.3-0.5 first ap- 
pear at the onset of convective instability around the radius 
of maximal neutrino heating, but before the dipolar mode of 
the shock reaches its limiting amplitude. These vortices are 
a source of acoustic waves, which have a similar period to 
the large-scale oscillation of the shock. Near the threshold 
for an explosion, the turbulence in the gain region becomes 
supersonic, as the existence of widespread secondary shocks 
attests. These shocks convert turbulent kinetic energy to in- 
ternal energy, increasing the effective heating rate. 

There are at least two reasons why explosions by the mech- 
anism investigated here may be more difficult in fully three- 
dimensional simulations. First, the existence of more degrees 
of freedom for the low-order modes of the shock in three di- 
mensions implies that the amplitude of individual shock oscil- 
lations is lower. As a result, it is more difficult for the shock to 
extend out to the radius where a-particle recombination gives 
it the final push. Second, an axisymmetric explosion that is 
driven by neutrino heating involves the accumulation of vor- 
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ticity on the largest spatial scales, an effect that is special to 
two dimensions. A full resolution of these issues is possible 
only with high-resolution three-dimensional simulations. 
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APPENDIX 



A. ALPHA-PARTICLE ABUNDANCE IN NUCLEAR STATISTICAL EQUILIBRIUM 

We calculate the a-particle mass fraction X a in nuclear statistical equilibrium by limiting the nuclear species to a-particles and 
free nucleons, and fixing the electron fraction Y e = 0.5. We tabulate X a and temperature T as a function of pressure p and de nsity 
p, and then use these tables to calculate the rate of release of nuclear binding energy by the method described in 32.1.11 The 
temperature does not appear explicitly in the FLASH hydrodynamic solver, and only enters the flow equations indirectly through 
X a . 

We include the contributions to p from radiation, relativistic and partially degenerate electron-positron pairs, and nonrelativistic 
a-particles and nucleons. When k B T > m e c 2 /2, it can be written (Bet he et"aii r980): 



P = 



1 (k B Tf 
12 (Jic) 3 



llTT 2 , 1 , 

—— + 2rf +— ^ 

15 7T / 



+ ( \ -^k B T, 
4 / m u 



(Al) 



where rj = p e /(k B T) the normalized electron chemical potential, also known as degeneracy parameter, and ft, c, and m u are 
Planck's constant, the speed of light, and the atomic mass unit, respectively. The density and degeneracy parameter are further 
related by 



where Y e is the electron fraction. The equilibrium fraction of a-particles is given by the nuclear Saha equation, 

/ Q a \ ( m u k B T} 3 >- 
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as supplemented by the conditions of mass and charge conservation, 

X n +X p +X a = 1 

1 

X„ + —X ri — Yf, . 



(A3) 

(A4) 
(A5) 



In eqs. (IA3t - (IA5t . X„ and X p are the mass fractions of free neutrons and protons, respectively, and Q a = 28.3 MeV is the binding 
energy of an a-particle. Combining eqs. ( IA1I ) and ( IA2b gives rj and T in terms of p and p. The equilibrium mass fraction X^ 
is calculated from p and T. For numerical calculations, we tabulate X^, dX ^ / di n p, a nd dX ^/dlnp for a grid of density and 
pressure. In addition, we tabulate partial derivatives of T to substitute in eqs. ( IB 10b and (IB 1 II) . 

Figure [22] shows contours of constant X^ and constant entropy for different varia bles as a function o f density. The entropy per 
nucleon is obtained by adding the contributions from the different components (e.g. Bet he et al.l [l980). 



rj(TT 2 +T] 2 ) 



■ + ln 



i u n Q {T) \ 
P J 



-X p \nXp-X n \nX n - -X a In (X a /32) . 



(A6) 



The postshock density in the initial configuration is typically p2 ~ 10 9 g cm" 3 with an entropy ~ 10- 15^/nucleon. The formation 
of a-particles that is seen in Figure [TJ results from an expansion of the shock into the part of the thermodynamic plane in 
Figs.l22h,b where pi < 10 9 g cm" 3 and T < 1 MeV. In this regime, the electrons are non-degenerate and the pressure in photons 
and pairs begins to exceed the nucleon pressure. The dip in the adiabatic index seen in Figure l22f results from a-particle 
dissociation/recombination, which partially compensates the change in internal energy due to compression/expansion. 



B. TIME-INDEPENDENT FLOW EQUATIONS DESCRIBING INITIAL MODELS 

We write down the ordinary differential equations that are used to compute the initial flow, and the density profiles in Fig- 
ure[2k,b,c. The steady state Euler equations in spherical symmetry are 



1 dv,. 1 dp 2 
-^ + -^ + -=0 
v r ar p ar r 

dv r 1 dp 

dr p dr 
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FIG. 22. — E quat i on of state of a fluid containing n, p, a, photons, and finite-temperature and partially degenerate electrons, in nuclear statistical equilibrium. 
We solve eqns. |AT)-(A5) and tabulate all quantities on a grid of density and pressure. Panels (a), (b) and (c): pressure, temperature, and degeneracy parameter as 
a function of density, for fixed values of and entropy. Dotted line: ri = n, the approximate boundary between degenerate and non-degenerate electrons. The 
gray area shows the region of the thermodynamic plane where X^ 1 > 0.5. Panel (d): partial derivatives of X^ 1 with respect to density (solid lines, positive) and 
pressure (dashed lines, negative). Panel (e): ratio of relativistic pressure (photons and pairs) to material pressure (a-particles and nucleons). Panel (f): adiabatic 
index 7 for different adiabats (dotted line: 7 = 4/3). 



dr p dr 

where e mt is the internal energy per unit mass, Jz?c, and S£ a the source terms described in eqns. d9l)-dTTb. and g = GM/r 2 . 
Since two variables suffice to describe the thermodynamic state of a system, we write 



deint _ E dp dp 
dr dr p dr 

and 

'dr +Ap dr 



Se a =A t -f. +A p -^. (B5) 



The c oeffi cient s Ej and A,- encode the dependence on the equation of state. Replacing eqns. (IB4t and (|B5t in (|B3t , and using 
eqns. ( IBU and (IB2l i to eliminate the pressure derivative, we obtain 

dp = (pv r Ep -A p ) [pg - 2pv 2 r /r) + (& H - S? c ) 
dr (pv r E p -pv r / p-A p )+v 2 r (pv r E p -A p ) 

The coefficients in eqs. (IB4b and ( |B5t work out to 

1 

( 7 -l)p 



E p = ^— (constant 7) (B7) 



e p = —, — ~7T^> (constant 7). (B8) 

For a constant-7 equation of state, e\ nl = p/[Qy- l)p]. The pressure (equation IA1II ) in the NSE model described in Appendix 
[A] can be decomposed into contributions from relativistic particles and from nucleons, p = p K \+p mat , and the specific internal 
energy is 

1/ 3 \ p 3 / 3 \k B T 

emt = - 3p re l + -Pmat = 3 ^ - - 1 - -X a . (B9) 

\ 2 / p 2 V 4 J m u 
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One therefore finds 

3 9 k B T dX a 3 / 3 \ 1 d(k B T) 

£ P = -+a— 1-jXa ^f 2 , (NSE) (BIO) 

F p 8 m u dp 2 \ 4 J m u dp 

3 P ^ kBTdXg 3 A_3^ \ 1 9(fc g r ) 

8 m„ dp 2 \ 4 Q / /«„ <9p 

The initial postshock solution is obtained by integrating the above equations from r s to an inner radius r» at which the flow 
stagnates. We iterate the normalization of the cooling function in equation (0 so that r* = 0.4r s o in the absence of heating. When 
adding heating, the cooling normalization and r* are kept fixed, which results in an expansion of the shock from its initial position 
to r s > r s o (Figure |2^,b). The initial Mach number at the inner boundary is chosen so as to satisfy 



E p = — I+c^^T-ol 1 "!^)— J ^T L - ( NSE ) < B11 ) 



(-2H.; _ -Sfc,! + -Sfa I /) Vi 



: 0.995 



\M\, (B12) 



where the sum is taken over the computational cells below the shock at our fixed resolution (see j32.1.3l . V,- is the volume of 
each computational cell, and the source terms are evaluated at the inner radial cell face. The numerical coefficient on the right 
hand side depends on the radial resolution, and is chosen empirically to prevent runaway cooling due to discreteness effects in 
time-dependent calculations. The resulting inner Mach number is typically 10~ 3 - 10~ 2 . 

When including a-particles in the EOS, one needs to calculate self-consistently the value of below the shock, the corre- 
sponding dissociation energy e(t = 0) [equation <j8j] , and compression factor k [equation Q]. The density upstream of the shock 
is obtained from 

P^= a 2 M 1 m » (B13) 
Airri\v\{r s )\ 

where 

n(r s ) = (B14) 

<l+2Mf(r s )/h-l) 



is the upstream velocity at r = r s , while the upstream pressure satisfies 

P i(^)[vi(^) ] 2 
lM 2 Ar s ) 



Eqns. ( IB 13b and ( IB 15b are transformed to physical units for input to the NSE model by adopting Afi(r s o) = 5, M = Q3Mq s 
M = 1 .3M Q , and a particular value for the shock radius r j0 in the absence of heating. 

C. TIME EVOLUTION 

Here we give some further details of the time evolution of our initial models using FLASH2.5. Heating and cooling are applied 
in an operator split way in between hydrodynamic sweeps. Nuclear dissociation in the constant-e model is implemented through 
the fuel+ash module in FLASH, with the modifications described in Paper I. The rate of change of specific internal energy is 
computed using the current hydrodynamic variables and timestep, after which the EOS subroutine is called to ensure that the 
variables are thermodynamic ally consistent. 

In the NSE nuclear dissociation module, numerical stability is maintained using an implicit update of the pressure in between 
hydro sweeps, 

Pnew = Pcm + (j- l)P<? nuc (p, Aiew), (CI) 

where e nuc is the energy generation per timestep in equation ©, and the subscripts cur and new refer to the current and new 
value of the pressure, respectively. The density is kept constant across this step, so as to be consistent with the other source 
terms. Equation dC It usually converges in 3-4 Newton iterations, adding a negligible overhead to our execution time. We 
restrict the timestep of the simulation so that, in addition to the standard Courant-Friedrichs-Levy condition, it enforces |e nuc | < 
0.8(jt>//5)/(7— 1). To prevent a-particle recombination in the cooling layer (due to the decrease of internal energy), we adopt a 
cutoff in density, so that X a = if p > 3 x 10 10 g cm -3 . 

Both the constant-e and the NSE dissociation modules require that the Mach number remain below a fixed value Alburn for 
burning, which has the effect of preventing dissociation or recombination upstream of the shock (Paper I). This results in a small 
amount of incomplete burning in the presence of strong shock deformations, a phenomenon which is also encountered in the 
full collapse problem. We set the threshold to Alburn = 2 in most of our simulations. The critical heating parameter H a depends 
weakly on A^bum : changes in Album cause small changes in the amount of unburnt material with zero Bernoulli parameter, and 
only slightly alters the net energy of the gain region. At the outer boundary of the simulation volume, AJburn is just below the 
Mach number of the upstream flow. We have tried expanded the outer boundary to r = 9r S Q (with a somewhat smaller Album) an d 
found that runs that did hit r = 7r s o still hit the new outer boundary. 

Aside from the inclusion of heating and NSE dissociation, the numerical setup for our runs is identical to that of Paper I. In 
that paper, the numerical output was verified by comparing the measured growth rates of linearly unstable modes of the shock 
with the solution to the eigenvalue problem. We have tested the implementation of our NSE dissociation model by verifying that, 
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in the absence of initial perturbations, our steady state initial conditions remain steady. Spherical transients present in the initial 
data die out in a few I = oscillation cycles, and are present even when nuclear burning is omitted. (See Paper I for a more 
extended description.) 
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